function IdSet=OneSide3_ind(u01grid, u11grid, u21grid,u02grid, u12grid, u22grid,...
                            sgb,equalorder, ...
                            Ugridreduced, sgreduced,...
                            id_Ugridreduced,...
                            s_id_gridreduced, sU_id_gridreduced,...
                            n_U_sets, Tcoord, indices_pairs,...
                            numeqconstraints,coordrowseq, ...
                            d3, e5, f5, g3, h5, i5, l3, m5, n5, ...
                            fillAeq, beq)
%% Loop over parameter values 
check_exists=zeros(sgreduced,1);
eps=10^(-6);
param_MOSEK.MSK_IPAR_LOG = 0; 
param_MOSEK.MSK_IPAR_INTPNT_BASIS = 'MSK_BI_NEVER';
param_MOSEK.MSK_IPAR_NUM_THREADS = 1;
for g=1:sgreduced
% Set indices to keep track of equal elements     
idU1=id_Ugridreduced{1}(g,:);
idU2=id_Ugridreduced{2}(g,:);
idU3=id_Ugridreduced{3}(g,:);
idU1_t=id_Ugridreduced{1}(g,:).';
idU2_t=id_Ugridreduced{2}(g,:).';
idU3_t=id_Ugridreduced{3}(g,:).';
% Relevant sizes of the vector of unknowns
s1=s_id_gridreduced(g,1);
s2=s_id_gridreduced(g,2);
%s3=s_id_gridreduced(g,3);
sU=sU_id_gridreduced(g); 



%%%%%%%%%%%%%%%%%%%%%% Equality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Coordcolumnseq
id_01=[changecoord(s1,s2,idU1((1:15+3*2*sgb)),idU2(15),idU3(1))...  
       changecoord(s1,s2,idU1(15),idU2((1:14)),idU3(1))...
       changecoord(s1,s2,idU1(15),idU2(16:15+3*2*sgb),idU3(1))...
       changecoord(s1,s2,idU1(repmat((1:1:15+3*2*sgb),1,13+2*sgb*3)), idU2(15), idU3(kron((3:1:15+3*2*sgb), ones(1,15+3*2*sgb))))...
       changecoord(s1,s2,idU1(15),idU2(repmat((1:14),1,13+2*sgb*3)),idU3(kron((3:1:15+3*2*sgb), ones(1,14))))...
       changecoord(s1,s2,idU1(15),idU2(repmat((16:15+3*2*sgb),1,13+2*sgb*3)),idU3(kron((3:1:15+3*2*sgb), ones(1,2*sgb*3))))...
       changecoord(s1,s2,idU1(repmat((1:1:15+3*2*sgb),1,15+3*2*sgb)),idU2(kron((1:1:15+3*2*sgb), ones(1,15+3*2*sgb))),idU3(2))...
       changecoord(s1,s2,idU1(14),idU2(14),idU3(1))]; %zeros and ones (row)


coordcolumnseq= [changecoord(s1,s2,idU1(1),idU2(14),idU3(3)) changecoord(s1,s2,idU1(14),idU2(14),idU3(3)) changecoord(s1,s2,idU1(1),idU2(14),idU3(1)) ... %observational equivalence
                 changecoord(s1,s2,idU1(14),idU2(14),idU3(3)) changecoord(s1,s2,idU1(14),idU2(1),idU3(3))...
                 changecoord(s1,s2,idU1(1),idU2(1),idU3(1)) ...
                 changecoord(s1,s2,idU1(7),idU2(14),idU3(9)) changecoord(s1,s2,idU1(14),idU2(14),idU3(9)) changecoord(s1,s2,idU1(7),idU2(14),idU3(1)) ... 
                 changecoord(s1,s2,idU1(14),idU2(14),idU3(9)) changecoord(s1,s2,idU1(14),idU2(7),idU3(9))...
                 changecoord(s1,s2,idU1(7),idU2(7),idU3(1)) ...
                 id_01 ... %zeros and ones
                 changecoord(s1,s2,idU1(1),idU2(14),idU3(1)) changecoord(s1,s2,idU1(2),idU2(14),idU3(1)) ... %symmetry
                 changecoord(s1,s2,idU1(3),idU2(14),idU3(1)) changecoord(s1,s2,idU1(4),idU2(14),idU3(1)) ... 
                 changecoord(s1,s2,idU1(5),idU2(14),idU3(1)) changecoord(s1,s2,idU1(6),idU2(14),idU3(1)) ... 
                 changecoord(s1,s2,idU1(7),idU2(14),idU3(1)) changecoord(s1,s2,idU1(8),idU2(14),idU3(1)) ... 
                 changecoord(s1,s2,idU1(9),idU2(14),idU3(1)) changecoord(s1,s2,idU1(10),idU2(14),idU3(1)) ... 
                 changecoord(s1,s2,idU1(11),idU2(14),idU3(1)) changecoord(s1,s2,idU1(12),idU2(14),idU3(1)) ... 
                 changecoord(s1,s2,idU1((16:1:15+3*2*sgb)),idU2(14),idU3(1))...
                 changecoord(s1,s2,idU1(14),idU2(1),idU3(1)) changecoord(s1,s2,idU1(14),idU2(2),idU3(1)) ... 
                 changecoord(s1,s2,idU1(14),idU2(3),idU3(1)) changecoord(s1,s2,idU1(14),idU2(4),idU3(1)) ... 
                 changecoord(s1,s2,idU1(14),idU2(5),idU3(1)) changecoord(s1,s2,idU1(14),idU2(6),idU3(1)) ... 
                 changecoord(s1,s2,idU1(14),idU2(7),idU3(1)) changecoord(s1,s2,idU1(14),idU2(8),idU3(1)) ... 
                 changecoord(s1,s2,idU1(14),idU2(9),idU3(1)) changecoord(s1,s2,idU1(14),idU2(10),idU3(1)) ... 
                 changecoord(s1,s2,idU1(14),idU2(11),idU3(1)) changecoord(s1,s2,idU1(14),idU2(12),idU3(1)) ... 
                 changecoord(s1,s2,idU1(14),idU2((16:1:15+3*2*sgb)),idU3(1))...   
                 changecoord(s1,s2,idU1(14),idU2(14),idU3(3)) changecoord(s1,s2,idU1(14),idU2(14),idU3(4)) ... 
                 changecoord(s1,s2,idU1(14),idU2(14),idU3(5)) changecoord(s1,s2,idU1(14),idU2(14),idU3(6)) ... 
                 changecoord(s1,s2,idU1(14),idU2(14),idU3(7)) changecoord(s1,s2,idU1(14),idU2(14),idU3(8)) ... 
                 changecoord(s1,s2,idU1(14),idU2(14),idU3(9)) changecoord(s1,s2,idU1(14),idU2(14),idU3(10)) ... 
                 changecoord(s1,s2,idU1(14),idU2(14),idU3(11)) changecoord(s1,s2,idU1(14),idU2(14),idU3(12)) ... 
                 changecoord(s1,s2,idU1(14),idU2(14),idU3(13)) changecoord(s1,s2,idU1(14),idU2(14),idU3(14)) ... 
                 changecoord(s1,s2,idU1(14),idU2(14),idU3((16:1:15+3*2*sgb)))...   
                 changecoord(s1,s2,idU1(13),idU2(14),idU3(1)) ...
                 changecoord(s1,s2,idU1(14),idU2(13),idU3(1)) ...
                 changecoord(s1,s2,idU1(14),idU2(14),idU3(15)) ...   
                 changecoord(s1,s2,idU1(1),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(3),idU3(1)) ... %identical marginals
                 changecoord(s1,s2,idU1(1),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(5)) ... 
                 changecoord(s1,s2,idU1(2),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(4),idU3(1)) ... 
                 changecoord(s1,s2,idU1(2),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(6)) ... 
                 changecoord(s1,s2,idU1(3),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(1),idU3(1)) ... 
                 changecoord(s1,s2,idU1(3),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(7)) ... 
                 changecoord(s1,s2,idU1(4),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(2),idU3(1) )... 
                 changecoord(s1,s2,idU1(4),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(8)) ... 
                 changecoord(s1,s2,idU1(5),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(5),idU3(1)) ... 
                 changecoord(s1,s2,idU1(5),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(3)) ... 
                 changecoord(s1,s2,idU1(6),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(6),idU3(1)) ... 
                 changecoord(s1,s2,idU1(6),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(4)) ... 
                 changecoord(s1,s2,idU1(7),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(9),idU3(1)) ... 
                 changecoord(s1,s2,idU1(7),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(11)) ... 
                 changecoord(s1,s2,idU1(8),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(10),idU3(1)) ... 
                 changecoord(s1,s2,idU1(8),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(12)) ... 
                 changecoord(s1,s2,idU1(9),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(7),idU3(1)) ... 
                 changecoord(s1,s2,idU1(9),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(13)) ... 
                 changecoord(s1,s2,idU1(10),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(8),idU3(1)) ... 
                 changecoord(s1,s2,idU1(10),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(14)) ... 
                 changecoord(s1,s2,idU1(11),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(11),idU3(1)) ... 
                 changecoord(s1,s2,idU1(11),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(9)) ... 
                 changecoord(s1,s2,idU1(12),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(12),idU3(1) )... 
                 changecoord(s1,s2,idU1(12),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(10)) ... 
                 changecoord(s1,s2,idU1(13),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(13),idU3(1)) ...
                 changecoord(s1,s2,idU1(13),idU2(14),idU3(1)) changecoord(s1,s2,idU1(14),idU2(14),idU3(15)) ...
                 changecoord(s1,s2,idU1(d3(:).'), idU2(e5(:).'), idU3(f5(:).')) ...
                 changecoord(s1,s2,idU1(g3(:).'), idU2(h5(:).'), idU3(i5(:).')) ...
                 changecoord(s1,s2,idU1(l3(:).'), idU2(m5(:).'), idU3(n5(:).'))];
             
Aeq=sparse(coordrowseq, coordcolumnseq, fillAeq, numeqconstraints, sU); 

%%%%%%%%%%%%%%%%%%%%%%  Inequality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Repeat the way we constructed Tcoord_specific but now using the parameter values under consideration
vectors = cellfun(@(x) {x(g,:)}, Ugridreduced); %cell 1xn_U_sets
%For j=1,...,n_U_sets, it picks the g-th row of Ugridreduced{j}
T_temp = cell(1,n_U_sets);
[T_temp{:}] = ndgrid(vectors{:}); 
T_temp = cat(n_U_sets+1, T_temp{:});
T = reshape(T_temp,[],n_U_sets); %matrix of dimension as Tcoord reporting the elements of all the possible n_U_sets-tuples from the U-sets

D=[T(indices_pairs(:,1),:) T(indices_pairs(:,2),:) ...
   Tcoord(indices_pairs(:,1),:) Tcoord(indices_pairs(:,2),:)]; %n_U_sets-tuples n_U_sets-tuples coordinates coordinates
    
%Save in D1 the rows of D_temp where the first triplet is <= then the second
%Only the columns 7:end (reporting coordinates) are relevant
D1=D(sum(D(:,1:3)<=D(:,4:6),2)==3, 7:end);
sd1=size(D1,1);
%Remember <= means: (<,<,<) or (<,=,=) or (<,<,=) or (=,=,=) [any order]

%Save in D2 the rows of D_temp where the first triplet is >= (excluding = = =) then the second
%Only the columns 7:end (reporting coordinates) are relevant
D2=D((sum(D(:,1:3)>D(:,4:6),2)==3) |... %>,>,>
       (sum(D(:,1:3)>D(:,4:6),2)==2 & sum(D(:,1:3)==D(:,4:6),2)==1) | ... %>,>,=
       (sum(D(:,1:3)>D(:,4:6),2)==1 & sum(D(:,1:3)==D(:,4:6),2)==2), 7:end); %>,=,=
sd2=size(D2,1);

coordrowsineq=[kron((1:1:sd1), ones(1,8))  kron((sd1+1:1:sd1+sd2), ones(1,8))]; %1x[(sd1+sd2)*8]
%8 comes from the number of all possible vertices generated by each pair of triplets

coordcolumns1=[changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,2)),idU3_t(D1(:,3))) ... 
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,2)),idU3_t(D1(:,3))) ...
                   changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,5)),idU3_t(D1(:,3))) ...
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,5)),idU3_t(D1(:,3))) ...
                   changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,2)),idU3_t(D1(:,6))) ...
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,2)),idU3_t(D1(:,6))) ...
                   changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,5)),idU3_t(D1(:,6))) ...
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,5)),idU3_t(D1(:,6)))]; %sd1x8

coordcolumns2=[changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,2)),idU3_t(D2(:,3))) ... 
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,2)),idU3_t(D2(:,3))) ...
                   changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,5)),idU3_t(D2(:,3))) ...
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,5)),idU3_t(D2(:,3))) ...
                   changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,2)),idU3_t(D2(:,6))) ...
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,2)),idU3_t(D2(:,6))) ...
                   changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,5)),idU3_t(D2(:,6))) ...
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,5)),idU3_t(D2(:,6)))]; %sd2x8
          
coordcolumnsineq=[reshape(coordcolumns1.', 1, 8*sd1) reshape(coordcolumns2.', 1, 8*sd2)];   %transform coordcolumns1,2 in a row vector by reshaping row-wise  
%1x[(sd1+sd2)*8]


fillAineq=[repmat([1 -1 -1 1 -1 1 1 -1], 1, sd1) repmat([-1 1 1 -1 1 -1 -1 1], 1, sd2)]; %remember: signs are flipped because inequalities are <= 
%1x[(sd1+sd2)*8]

bineq=zeros(sd1+sd2, 1);

Aineq=sparse(coordrowsineq, coordcolumnsineq, fillAineq, sd1+sd2, sU);   %[#inequalityconstraints]x[sU]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Run LP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% MOSEK
lower=zeros(sU,1)+eps;
upper=ones(sU,1)-eps;
lower(id_01)=0;
upper(id_01)=1;
prob.blx=lower; %lower bound unknowns
prob.ulx=upper; %upper bound unknowns
prob.c=zeros(sU,1); %objective function
prob.a=[Aeq; Aineq];
prob.blc=[beq; -Inf(size(Aineq,1),1)];
prob.buc=[beq; bineq];
[~,result]     = mosekopt('minimize echo(0)',prob, param_MOSEK);
if isfield(result,'sol')
if strcmp(result.sol.itr.prosta, 'PRIMAL_AND_DUAL_FEASIBLE') %if the primal is feasible
   check_exists(g)=1; 
end
end
end

%Enlarge results to the original grid of parameter vectors
check_exists=check_exists(equalorder);
IdSet = [u01grid(logical(check_exists),:)  u11grid(logical(check_exists),:)  u21grid(logical(check_exists),:) ...
         u02grid(logical(check_exists),:)  u12grid(logical(check_exists),:)  u22grid(logical(check_exists),:)];

